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The structured-population model has been widely used to study the spatial transmission of epidemics in 
human society. Many seminal works have demonstrated the impact of human mobility on the epidemic 
threshold, assuming that the contact pattern of individuals is mixing homogeneously. Inspired by the recent 
evidence of location-related factors in reality, we introduce two categories of location-specific 
heterogeneous human contact patterns into a phenomenological model based on the commuting and 
contagion processes, which significantly decrease the epidemic threshold and thus favor the outbreak of 
diseases. In more detail, we find that a monotonic mode presents for the variance of disease prevalence in 
dependence on the contact rates under the destination-driven contact scenario; while under the 
origin-driven scenario, enhancing the contact rate counterintuitively weakens the disease prevalence in 
some parametric regimes. The inclusion of heterogeneity of human contacts is expected to provide valuable 
support to public health implications. 

To study the epidemic spreading of viruses and diseases in human society, various mathematical models have 
been proposed during the past decades 13 . In particular, the network-based models have attracted a great deal 
of attention from diverse disciplines 4 9 . Plenty of works 1016 explored the transmission of etiological agent on 
networks, where nodes are individuals and links describe the contacts between individuals to transmit the 
infection. 

Using the particle-network framework to mould the entire system into structured populations, we arrive at an 
important class of models in modern epidemiology, namely, metapopulation. In this model, individuals live in 
discrete subpopulations and may transfer between subpopulations via the mobility pathways (or links) 17 26 . The 
reaction-diffusion or reaction- commuting mechanism is harnessed to sketch human daily contact and mobility. 
More precisely, the disease prevails inside each subpopulation (with homogeneously mixed as assumed), and 
transmits between subpopulations through the travel of infected individuals. Applied in simulating the spatial 
transmission at large geographic scale, this model manifests its power in predicting the pandemic outbreak and 
evaluating the effectiveness of intervention strategies 27 34 . 

Many works 17,20,21,23-24 have focused on the impact of human mobility on the threshold of disease outbreak, 
which generalized the reaction-diffusion process to deal with heterogeneous networking populations, and 
assumed the mobility scheme as a Markovian memoryless migration. Recent empirical findings on human 
mobility 35-3 * report that the commuting behaviors, characterized by individual recurrent movements between 
connected locations, dominates the human daily transportation. Refs. 25, 26 extended the metapopulation 
framework by introducing the element of recurrent commuting, which assumes that individuals have the memory 
of their original resident subpopulation and displaced commuters who stay at the 'destination' subpopulation 
cannot diffuse to other places but return to their residence with a certain rate. In contrast to the random diffusion 
scenario, the commuting system exhibits a novel epidemic threshold as well as the phenomenon of the saturation 
of propagation velocity. 

It is well known that the contact pattern of individuals dramatically impacts the spatiotemporal dynamics of 
epidemic spreading in a population 13,40 . The evolution of the epidemic process is characterized by the force of 
infection, which identifies the probability of acquiring infection for a susceptible individual due to his contacts 
with infectious ones. The force of infection is defined proportional to the density of infectious individuals, the 
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infection probability via a contact between susceptible (S) and infec- 
tious (I) individuals, and the contact pattern between individuals 1 . 
Since human contact pattern is often reflected in the disease trans- 
mission rate, which is usually assumed constant 1,41,42 , in metapopula- 
tion models, the general assumption that the disease transmission 
rate in all subpopulations is identical indicates the fact that the indi- 
vidual contacts inside each subpopulation are uniform. 

Within each subpopulation, the assumption of homogeneous con- 
tact (homogeneously mixing) is supported with recent empirical 
findings. For instance, digital equipments such as Bluetooth 43 , active 
Radio Frequency Identification (RFID) devices 44 " 46 , wireless sen- 
sors 47 , and WiFi 48,49 , have been applied to collect the data of human 
close proximity contacts in various social environments, which 
unveils that the distribution of the number of distinct persons that 
each individual is in contact with per day has a small squared coef- 
ficient of variation (CV 2 ) 45 ~ 49 . Therefore, human contacts inside each 
subpopulation may possess the characteristic contact rate reflecting 
the number of distinct persons each individual encounters per day. 

At different subpopulations, human contact pattern may present 
evident discrepancies. The location-specific factor has been reported 
as a potential driver inducing the substantial variation of disease 
incidence between populations in reality 50,51 . Generally, because of 
the distinction in social conditions or lifestyles, the contact rate of 
individuals living in an urban area is largely different from what we 
can expect for small-town residents. Therefore, we consider the vari- 
ety of contacts between connected subpopulations to study their 
impact, and in this paper, we theoretically analyze this issue with 
the phenomenological model proposed in Ref. 52, where the recur- 
rent commuting of individuals couples two typical subpopulations. 
The standard susceptible-infectious-susceptible (SIS) model is con- 
sidered to reflect individuals' compartment transitions. 

We first consider a general case that the individual contact features 
depend on the location where an individual is. This scenario illus- 
trates the influence of social environment of the located subpopu- 
lation, and may reflect human adaptive ability 53 55 . Intuitively, con- 
sidering the commuting flows between the urban and the suburb 
areas, one can expect that the contact rate for an individual commut- 
ing from the suburb to the urban area will increase, and vice versa. 
Therefore, we assign each individual in the located subpopulation x 
with the same characteristic contact rate c x , which means that on 
average he contacts with c x other individuals in the same subpopula- 
tion per unit time. We define this destination- driven mode as the 
type-I contact scenario. 

As the commuting mobility distinguishes individuals' original 
residences from their destinations, we further consider another case 
that the contacts of individuals correlate to their own residences. This 
scenario may derive from the anchoring effect of human 56,57 , and 
reflects that many social factors of one's home location, e.g., eco- 
nomy, cultural background, inevitably affect how gregarious a per- 
son is. For instance, the schools of related disciplines of a university 
are usually located in the same campus (here each campus is a sub- 
population), while each day students from different campuses might 
commute among the campuses by the school buses. Students 
majored in humanities or social sciences usually have a tendency 
to participate in social activities, and they may have a high contact 
rate; whereas students majored in natural sciences or medicine (par- 
ticularly, graduate students) are more prone to spend their time in 
the laboratories or libraries, and they may own a low contact rate. 
Despite students from different disciplines having distinct contact 
rates may reside at different campuses, the commuting movements 
still induce their intersectional couplings in each campus per day. For 
simplicity, we assume that each individual from subpopulation x has 
the characteristic contact rate c' x . This means that each one from x 
will on average have contacts with c' x other individuals per unit time, 
no matter where he locate at that time. We define this origin-driven 
mode as the type-II contact scenario. 



In both two cases, the force of infection for a susceptible stems 
from all his contacts with infectious individuals per unit of time. By 
means of analytical arguments as well as extensive computer simula- 
tions, we demonstrate that these location-specific heterogeneous 
contact scenarios significantly decrease the epidemic threshold of 
the entire system, and thus favor the disease outbreaks in more broad 
parametric regimes. Under the destination-driven scenario, the vari- 
ance of disease prevalence displays a monotonic mode as the contact 
rates change; while the results are more unexpected under the origin- 
driven scenario: Enhancing the contact rate will weaken the disease 
prevalence in some parametric regimes. 

Results 

We first specify the mechanism of commuting mobility to transfer 
individuals between subpopulations. Consider two coupled subpo- 
pulations, x, y, each of which has a population size of the residents 
N* = N* = N. They act as the reaction places where the contagion 
process occurs. The model proceeds at discrete time steps, with the 
unit interval as one hour. Individuals leave their original resident 
subpopulation x(y) to visit the neighboring 'destination' subpopula- 
tion y(x) with a per capita diffusion rate o xy (a yx ) at each time step, 
and the displaced commuters will return to their residence x(y) with a 
per capita return rate x x (x y ) per unit time. For simplicity, we assume 
the detailed balance as: a xy = a yx = a, x x = z y = t, and Figure 1 
illustrates the commuting process. 

According to Refs. 26, 52, the residents of each subpopulation can 
be partitioned into two subgroups based on the locations where they 
actually are at time t (see Figure 1). Define N^itfiN^it)) the sub- 
group size of the individuals staying in their original resident sub- 
population x(y) at time f, and N xy (t)(N yx (t)) the subgroup size of the 
individuals from subpopulation x(y) presenting in the subpopulation 
y(x) at time r. Therefore, the population size of the residents of each 
subpopulation is N* = N xx ( f ) + N xy (t),N^ = N yy ( t) + N yx (f).We use 
the mean-field approximation to describe the commuting process 
with a set of rate equations governing the evolution of population 
dynamics 26,52 , which are discussed in Supplementary Information, 
and we have the following equilibria: 

N* N* 
N e i = ? N eq = 



l + o/x' xy 1 + t/o-' 



N* 

N e i = y - N etl -- 

yy i+ffA' yx 



n; 



(i) 



l+T/V 



with the condition a 1 x \ i.e., averagely, individuals will remain 
at their residence most of the time, which mimics that in reality only a 




Figure 1 | Schematic illustration of the commuting process between two 
subpopulations x, y. The cyan and orange arrows indicate the back and 
forth commuting flows. 
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small fraction of residents leave to travel. This condition facilitates 
the analysis by considering that per unit time each subpopulation 
x(y) has an effective number of residents N^(N^n interacting with 
the individuals in subgroup N yy (N xx ). 

The contagion process takes place inside each subpopulation, 
which is composed of the infection dynamics and the contact 
dynamics here. For the infection dynamics, we consider the standard 
SIS compartment model. At time r, each individual falls in one of the 
disease compartments: susceptible or infectious. The population size, 
e.g., N xx (t) {N yx (f)l can be divided into S^fXS^f)), f„(i)(^), 
which is the number of susceptible and infectious individuals in each 
subgroup, respectively. At each time step, if a susceptible individual 
contacts an infectious one, he may acquire the infection with trans- 
mission rate fS. An infectious individual recovers and becomes sus- 
ceptible again with rate v per unit time. 

We now specify the contact dynamics in more detail. In the case of 
the type-I (destination-driven) contact scenario, we assume that 
per unit time each individual staying in subpopulation x(y) has the 
same characteristic contact rate c x {c y ). Figure SI (a) in Supporting 
Information illustrates this case. In the unit interval, the infection 
probability, X, for a susceptible individual is determined by all his 
contacts with the infectious ones. Therefore, the force of infection for 
any susceptible in subgroup N xy or N yy at time t is 

;^(0=^(0=cAW=/^(^(0+^(0)/(n:;+n;;), (2) 

where A y (t) characterizes the probability that a random contact with 
any individual staying in subpopulation y will lead to the infection at 
time t. Similarly, the force of infection for any susceptible in sub- 
group N yx or N xx at time t is 

V« = = Mt) = ik x (i yx (t) +i„(t)) I (n%+n%) . (3) 

With the mean-field approximation, the contagion process in each 
subgroup under the type-I scenario can be delineated by a set of rate 
equations (see Eqs. (S5)-(S8) in Supplementary Information). 

We next extend to the type-II (origin-driven) contact scenario. Per 
unit time, an individual from subpopulation x(y) has the character- 
istic rate of contacts, c' x (^c'^j, which means that on average he 
encounters with c' x ( c 'y) persons, no matter where he is (see Figure 
Sl(b) for the illustration). With the mean-field approximation, the 
contagion process in each subgroup under the type-II scenario can 
also be described by the formalism of rate equations as Eqs. (S5)-(S8), 
except the expressions of the force of infection. The force of infection 
for any susceptible staying in subgroup N xy at time t obeys the fol- 
lowing relation 

c%{t) =pc' x (c'Mt)+c' y i yy (t)) /(<.n^+ c ;n;;), ( 4 ) 

where XJi) denotes the probability that a random contact with any 
individual staying in subpopulation y will lead to the infection. It is 
worth emphasizing that to normalize the fraction of contacts of the 
infectious ones in each subpopulation, the denominator of Eq.(4) 
should be the total number of contacts of the population instead of 
the total number of individuals staying in the subpopulation. 
Similarly, the force of infection for a susceptible presenting in other 
subgroups at time t are as follows: 

i' yy (t)=c^ y (t)=pc;[c%(t)+c' y i yy (t)) /( c ^+ c ;n;;), (5) 
= W*) = K { c 'Mt) + c' y Ut)) I (« + c> y N%) , (6) 

C(t)=c'X(t)=K(<Ut)+c' y I yx (t)) /{c'XA + C^i) . (7) 

Besides, we consider a baseline scenario that the contact rate in each 
subpopulation is equal to (c) = (c x + c y )/2, which corresponds to the 
traditional situation that the contact pattern in all subpopulations is 



uniform. In this case, the epidemic threshold is captured by the basic 
reproductive number, R' Q = fS{c) / v, which identifies the expected 
number of secondary infections produced by an infected individual 
during his infectious period within the entire susceptible popula- 
tion 1-3,41,58 . lfR' 0 <l, the whole system remains at the disease-free state. 

Analysis. With embedding the aforementioned contact features, the 
epidemic threshold can be characterized by an extended version of 
the basic reproductive number lZo 2A1 ' 5S - We first build the next 
generation matrix (NGM) TS -1 from the rate equations (S5)-(S8) 
at the disease-free equilibrium (DFE). Under the type-I contact 
scenario, the transmission matrix T is defined as 



(Br* 
dl„ 



T\Jxxi %yxf lyyi ^xy) 



8T XX \ 



di- 



xy 



ST, 



xy 



dT, 



xy 



dlxyj 



f {l-a)jic x x (l-g)feT 

T + ff T + O" 

(l — z)fic x cr (1 — T)fic x c 



(\—o)ficy% (l — a)PcyX 

T + O T + O" 

(1— x)pc y a (l—-z)pc y o 



while under the type-II contact scenario, the transmission matrix Tis 



(8) 



T {jxxi lyxi ^-yy> ^-xy) 



dlxx 



dT„ 



\SI XX 

( {\-a)pc' 2 x (l-g)fa' x c' ; 
c' x + c' y a/i: 

{\-X)Pc' x C'y 
C' X T/(7 + C'y 



dT xx \ 
8L, 



8T xy 
dl xy j 



C' X + C'yO/t 



/a + c 



\ 



(\~a)pd x c'y 
c' x o/x + c'y c' x a/-c + c' y 
(\-z)pc' x c'y (1-t)Pc' 2 x 
c' x + dy%j a c' x + c'yt/a J 



(9) 



In the above two matrices, each entry Tjj(i,j e (x, y)) represents the 
rate of generating new infections in subgroup Nu. In these two 
scenarios, the transition matrix 2 is 



dlv 



2(^ct> Iy X ? ^yy-> I X y) — 



y X , iyy, * X y ) 

( ff + (l-ff)v 
0 
0 



dl 



xy 



\ dl xx 
0 

t + (1-t)v 

— T 
0 



xy 



8Z X y 

~dlxy~ j 



0 

— a 

(l- ff )v 
0 



(10) 



0 
0 

(1-T>/ 
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where E,j = E ( ^~ — Et. Et denotes the rate of transferring infectious 
individuals into subgroup Njp while E^ denotes the rate of 
transferring infectious people out of Ny. The recovery of infective 
individuals also decreases the number of infectious ones, and thus is 
included in E^" . The matrices Tand E are obtained at the disease-free 
equilibrium. The global version of the basic reproductive number is 
calculated by the spectral radius of the NGM 41 ' 58 : 

fto^TE- 1 ). (11) 

Since the explicit expression of the threshold is unavailable, we 
present the numerical results in the following. 

Figure 2 presents the c x —c y (c' x —c'^j phase diagrams of the global 
1Z 0 . As the empirical evidence tells that the commuting takes place at 
the time scale t -1 ~ 1/3 day 36 , we choose the return rate x = 0.125/ 
hour. For the diffusion rate, we first set a large value a = 0.042/hour, 
which implies that the time scale of diffusion is 1 day. We set other 
parameters of the disease as: v = 0.042/hour (a short duration for the 
mean infectious period), and /i = 0.021/hour. Note that fixing the 
disease parameters as other values does not change the main results 
of this paper. The white dashed line in each panel indicates the 
parametric values corresponding to the threshold 1Z Q = 1. The evid- 
ent difference of the thresholds between the two heterogeneous con- 
tact scenarios shows that the type-II contact manner may induce the 
disease outbreak at much lower contact rates. This is due to the fact 
that the commuting process facilitates each subpopulation x(y) to 
share an effective number of individuals N xy {NyV) into the neighbor- 
ing subpopulation y(x). As the displaced commuters will closely 
interact with local residents, individuals with distinct contact rates 
have the opportunity to directly encounter in each subpopulation 
under the type-II scenario. 

As shown in Figure 2(a), one can find a monotonic mode for the 
variance of global IZo under the type-I scenario, which positively 
correlates to the increase of contact rates c x and c„. The results are 
more complex under the type-II scenario (see Figure 2(b)), because 
the values of global TZq at the upper-left and lower-right corners of 
the phase diagram are also very large. The variance of TZq does not 
show a monotonic mode when c' x or c' y is large (> 2.34) in the case of 
type-II scenario. For instance, when c' x — 3, TZg gradually decreases 
until reaches the bottom, then monotonically grows with the increase 

ofc;. 

We further compare the parametric regions of the endemic phases 
among the baseline, type-I, and type-II scenarios with the aforemen- 
tioned parameters, as shown in Figure 3(a). The upper triangular 
light orange area corresponds to the baseline with R' Q = /?(c) /v> 1. 
The parametric space of the type-I or type-II scenario with the global 
TZg > 1 is much broader than that of the baseline case. Therefore, 
the inclusion of both two location-specific heterogeneous contact 



patterns favors the epidemic outbreak. It is also of interest to inspect 
the effect of a smaller diffusion rate a on the threshold, since reducing 
a weakens the coupling between subpopulations. Figure 3(b) com- 
pares the parametric regions pertaining to the endemic phases 
among the baseline, type-I, and type-II scenarios with a = 0.021/ 
hour, which yields the same results. 

Simulations. To verify the above analysis, we use the dynamic Monte 
Carlo method to simulate the epidemic evolution. We focus on 
inspecting the parametric space of contact rates. With each set of 
parameters, we monitor the dynamical evolution of epidemic 
spreading in each subpopulation, and record the disease pre- 
valence (the fraction of infectious individuals) in subpopulation 
x(y), I x (t)/N x (I y (t)/N y ), where N x = N xx {t) + N yx {t) (N y = N„(t) 
+ N xy (t)). Due to the presence of stochasticity, each independent 
Monte Carlo simulation generates a random realization of the 
dynamic process. We perform 10 3 random Monte Carlo simula- 
tions for each set of parameters. 

Initially, each subpopulation has 10 5 individuals, i.e., N xx (0) = 
N yy (Q) = 10 5 , N xy (Q) = N yx (0) = 0. We fix the return rate % = 
0.125/hour and the diffusion rate a = 0.042/hour to stress the coup- 
ling effect. Figure S2 presents the time series of population evolution 
in each subgroup. As analyzed in Supplementary Information, the 
population size of each subgroup quickly reaches to the equilibrium: 
N xx =N yy 1 = 7.5x 10 4 , N e x ^=N y l = 2.5x 10 4 . Then we simulate the 
contagion process with a variety of typical sets of contact rates, fixing 
the transmission parameters \i = 0.021, v = 0.042/hour. The status of 
each individual is updated in parallel per unit time. The mobility and 
disease parameters are directly converted into the probabilities with 
the unitary time scale of one hour. The algorithm details of the 
simulation see the section of Methods. 

We start the contagion process by introducing the seeds of infec- 
tious individuals into a given subgroup. For simplicity, the case of 
multiple introductions is not considered here, i.e., the seeds are only 
introduced into one subgroup. The contagion process begins with 
I xx (0) = 50, Ste(0)=.Ng-jr„(0), I y M = I yy (0) = I xy (Q) = 0. By 
fixing the transmission parameters = 0.021, v = 0.042, we mainly 
study the stationary prevalence (SP) If /N°i, If /Nf with different 
parameters (SP reflects the average fraction of infectious individuals 
in each subpopulation at the equilibrium). To obtain each data point, 
we perform 10 3 Monte Carlo random experiments, each of which is 
simulated with 10 4 time steps. We define an outbreak run as the 
realization that the number of infectious individuals at the equilib- 
rium is larger than 50. The stationary prevalence is calculated by 
averaging over all the outbreak realizations. We first record the mean 
value of the number of infectious individuals in each subpopulation 
over the last 10% time steps of each outbreak simulation, then 
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Figure 3 | Comparison of the parametric space pertaining to the endemic phases among the baseline, type-I, and type-II contact scenarios, t = 0.125, 
v = 0.042, ft = 0.021/hour, and a = 0.042/hour in (a) and a = 0.021/hour in (b). 



gradually increases, fixing 



average them to obtain the stationary prevalence. If there does not the variance of I x q /N^, Iy q /Ny q as 1 

exist an outbreak run, I x q /N e x q = I eq /N e y q = 0. c x = c' x = 1. The thresholds of contact rate c y Ic'A , which separate the 

In Figs. 4-5, we compare the phase diagrams of the stationary transition from the disease-free phase to the endemic phase, of both 

prevalence f q j N e x q , f q j N e y q among the baseline, type-I, and type- the type-I and type-II scenarios are much smaller than that of the 

II contact scenarios with several typical c x (c' x ) . Figures 4(a)-(b) show baseline case, while the threshold of type-II scenario is the smallest. 
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Figure 4 | Comparison of the phase diagrams of j N^ q , ip/N^ among the baseline, type-I, and type-II scenarios. (a)(b) show the variance of 
I^/Nx 1 ^dI y q /Ny q as Cy\c y J gradually increases, fixing c x = c' x = 1 (c)(d) show the variance off x q /N x q and 7^ /Ny 1 as ^{^^j gradually increases, fixing 
c x = d x = 2. The error bar is not shown since the standard deviation is much smaller (less than 10~' of the mean value). The vertical colored dashed lines 
illustrate the theoretical prediction of the threshold of c y 
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Figure 5 | Comparison of the phase diagrams of I x q j N x q , I y q / N y q among the baseline, type-I, and type-II contact scenarios, (a) (b) show the variance of 
I? /N? andtf 1 as c y(^^J gradually increases, fixing c x = d x = 3. (c)(d) show the variance of f x /N x q and f q /N y q as c y ^cJ,J gradually increases, fixing 
c x = c' x = 4. The inset in (c) shows the enlargement of the yellow area in that panel. The error bar is not shown since the standard deviation is much smaller. 



We also observe the well agreement between the simulations and 
analysis, as illustrated by the colored dashed lines in each panel. In 
the endemic phase, we find a counterintuitive phenomenon that, 
with the same parameters, the stationary prevalence f x /N x q , 
Iy q /Ny q of the type-II scenario are larger than those of the type-I 
scenario. In subpopulation x, the larger f x q /N x q of the type-II scen- 
ario might relate to the introduction of individuals with a higher 
contact rate. However, in subpopulation y, it seems a little unexpec- 
ted that the type-II scenario leads to a larger f y q j Ny, because under 
the type-II scenario the displaced commuters from subpopulation x 
to y have a smaller contact rate when the global 7^o>l- In 
Figure 6(a), we compare the phase diagrams of global 72-o between 
the type-I and type-II contact scenarios, fixing c x = c' x =l. When 
7^0 > 1, with the same contact rate Cy = c' y , the global TZg of the 
type-II scenario is evidently larger than that of the type-I scenario, 
which means that with the same parameters, the prevalence is more 
serious in the former one. Since we assume that the individual con- 
tact rate keeps unchanged in the case of type-II scenario, this out- 
come may be due to the fact that individuals staying in subgroup N yy 
increase their "within self-subgroup" contacts compared with those 
of the type-I scenario. 

Figures 4(c)-(d) present the variance of stationary prevalence 
Ix q /N x q , f y q /N y q as Cy(c' y \ gradually increases, fixing c x = c' x = 2. 
In this case, the thresholds of contact rate c y (c'y^J f° r the baseline, 
type-I, and type-II scenarios overlap at c y — c' =2, which is also 



supported by Eq.(ll) as indicated by the red dashed line in each 
panel. In the endemic phase of 72-o > 1> the difference of stationary 
prevalence between the type-I and type-II processes is largely shrunk. 
Figure 6(b) compares the phase diagrams of global 7?.o with fixing 

c x — c' x — 2. In the case of a large c y \c' y j > the difference of 1Z Q becomes 
smaller than that of the former situation c x = c' x = \. 

When we further increase c x (c' x ), e.g., c x = c^. = 3 or c x = c' K = 4, 
the outcomes are more complex. Figures 5(a)-(b) show the variance 

of f q j N x q , Iy q j N y q as Cy [c'yj gradually increases, fixing c x = c' x = 3. 
In the baseline case, both the simulations and analysis agree that the 
threshold of contact rate decreases to c y = 1. The heterogeneous 
contact scenarios yield that the threshold of Cy\c'y\ vanishes. 

Figures 5(c)-(d) present the variance of ff /N x q , f y q /N e y q in depend- 
ence on c y (c'y^j > fixing c x = c' x = 4. Since the threshold of contact rate 

of the baseline also vanishes, the baseline results are not included 
in Figures 5(c)(d). The analytical results of global IZo (see 
Figures 6(c) (d)) indicate that the disease-free phase is eliminated, 
no matter how small c y ( c 'y) i s - 

Moreover, under the type-I process, both f q /Nx and I y q /N^ 
monotonically increase with the augment of c y , as shown in 
Figure 5. Under the type-II process, it is evident that I y q /Ny q also 
monotonically increases with the growth of c', while f q /N x q 
first gradually decreases, after reaching the bottom it begins to 
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Figure 6 | Comparison of the phase diagrams of global IZa between the type-I and type-II contact scenarios. Each panel shows the variance of global IZo 
as Cy{d^j gradually increases, with c x = d x = 1 in (a), c x = d x = 2 in (b), c x = d x = 3 in (c), c x = d x = 4 in (d). The horizontal gray dashed lines in the upper two 
panels illustrate the threshold value 7?.o = 1 • 



monotonically grow with the increase of c' y . The variance of global 
T^o under the type-II scenario also experiences a valley as c' y gradually 
increases (see Figure 6). To explain this untrivial phenomenon, we 
observe the stationary prevalence of each subgroup in detail. As an 
example, we present the results of c' x — 3 in Figure 7. It is clear that 
Iyy j Nyy and Iy X j Ny X monotonically increase with the increment of 
c' y , whereas the variances of f x \ j N xx and I X y j N X y are nonmonotonic. 
In subpopulation x, the probability of an individual of subgroup N xx 
to involve in a contact is c' x j [c' x N xx + c' y Ny X ^j , while the probability 
of an individual of subgroup N yx to involve in a contact is 
c 'y j ( c ^ xx c 'y^y*) • 'With the increase of c' y , the individuals in 
subgroup N xx will reduce their 'within self-subgroup' contacts. As 



c' y initially grows from 0.1, the net loss of infectious individuals in 
subgroup N xx suppresses the gain obtained from subgroup N yx , 
because N^^N^. In this stage, the presence of commuters from 
subpopulation y mitigates the epidemic situation in subpopulation x, 
and their influence will be enhanced with the increment of c' y . Thus 
one can observe the counterintuitive drop of f x /N?, In subpopula- 
tion y, with the increase of c' y , the net increase of infectious indivi- 
duals in subgroup N yy can compensate the loss of subgroup N xy , 

because N^>N^. Therefore, f y q /Ny q keeps on increasing as shown 
in Figure 5(b). As subpopulations x, y are coupled together, increas- 
ing c' y will eventually cease the decline of ff/N^. After hitting the 
bottom, f x j~bf x positively correlates with the increase of c' y . 
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Figure 7 | Variance of the stationary prevalence of infectious individuals in subgroups iV„, N v (a), JV^., Nyx (b) as c y y:' y j gradually increases, under 
the type-II scenario. We fix d = 3. 
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Figure 8 | The c x —c y yc' x —c' y j phase diagrams of I x q /N x q , I y q /N y q . (a)(b) show the type-I scenario. (c)(d) show the type-II scenario. The white dashed 
curve in each panel illustrates the threshold of contact rates leading to the global IZo = 1 . 



As a brief summary, we present the c x —c y (c' x —c'^j phase dia- 
grams of the stationary prevalence f x /N x q and f y q /N y q under the 
type-I and type-II scenarios, and give a holistic view about the impact 
of different contact patterns on the spatial transmission between 
populations. Figures 8(a)-(b) plot the c x —c y phase diagrams of 
I x q /N x q and f q /N y q under the type-I scenario, respectively. One 
can clearly observe the monotonic increase of f x /tf q and f q /N y q 

with the increment of c x and c y . The white dashed line in each panel 
illustrates the threshold of contact rates leading to the global IZo — 1 
(calculated by the NGM method), which successfully separates the 
endemic phase from the disease-free phase. Similarly, Figures 8(c)- 
(d) show the c ! x —c' y phase diagrams of I x q /N x q and f q / N y q under the 
type-II scenario, respectively. The parametric space corresponding to 
the endemic phase under the type-II scenario is broader than that of 
the type-I scenario. When c' x or c' y >2.34, the stationary prevalence 
f q I N e x q and f q j N y q does not always monotonically increase with the 
augment of c' x , c' y . In this case, the increment of contact rates will 
reduce the disease prevalence f q j N x q or I y q /Ny in some parametric 
regimes. 

Discussion 

With the evidence of location-related factors in reality, we have 
introduced two categories of location-specific contact patterns in a 
phenomenological structured populations model based on the com- 
muting and contagion processes. Through theoretical analysis and 
extensive computational simulations, we have shown that these het- 
erogeneous contact scenarios favor the disease outbreaks and evi- 
dently decrease the epidemic threshold of the entire system. This 
finding is robust against the changing of system parameters. More 
specifically, under the destination-driven scenario, the variance of 
disease prevalence exhibits the monotonic dependence of contact 
rates, while under the origin-driven scenario the increment of 



contact rates unexpectedly reduces the disease prevalence in some 
parametric regimes. Compared to the traditional framework, where 
the contact pattern in all subpopulations is uniform, the models 
presented in this paper demonstrate that the inclusion of these extra 
factors raises the risk of infection, mainly due to the presence of 
couplings between populations. Therefore, the interventions to limit 
the mobility flows, e.g., entry screening, and travel restriction, may 
benefit to control the disease invasion. To better understand the 
impact of human contact behaviors on the spatial transmission of 
infectious disease, more efforts are deserved to deal with more het- 
erogeneous networking populations in further studies. 

Methods 

Algorithm details. The contagion and commuting processes at each time step 
proceed as follows: 

(i) Contagion process. For the type-I scenario, the probability of a susceptible 
staying in subpopulation x to acquire the infection by each of his contacts at time t is 
A x (t) = jSCWO + + N yx {t)). Therefore, the force of infection for each 

susceptible staying in subgroup N** or N yx is ^(f) = A yx (t) = c x l x (t). Similarly, the 
force of infection for each susceptible staying in subgroup N yy or N xy at time t is A yy {t) 
= A xy (t) = c y A y (t). At this time step, the number of new infections in subgroup N xy is 
generated from a binomial distribution with probability A xy (t) and the number of 
trials S xy (t). The number of new infections in other subgroups is generated in the same 
way. At the same time, any infectious individual recovers and becomes susceptible 
again with rate v. The number of recovered individuals in each subgroup is generated 
from a binomial distribution with probability v and the number of trials defined by 
the number of infectious individuals in that subgroup. 

For the type-II scenario, each resident of subpopulation x(y) contacts with c' x ( c' J 

other individuals in his current location per unit time. At time f, the probability of a 
susceptible presenting in subpopulation y to acquire the infection by each of his 

contacts is X' y {t) = p(c' x I xy (i)+c' y I ry (t)\ j [c' x N xy (t) + c' y N„(t)^ . Thus the force of 
infection of each susceptible in subgroup Nxy is ^™(t) = c' x A y (t), while the force of 
infection of each susceptible in subgroup is ^{t) = c' y A.' y (t). The forces of 
infection of susceptible individuals in other subgroups are obtained similarly. At this 
time step, the number of new infections in subgroup N xy is generated from a binomial 
distribution with probability A' (f) and the number of trials 5^(t). The number of 
new infections in other subgroups is generated in the same way. At each time step, the 
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number of recovered individuals in each subgroup is generated from a binomial 
distribution with probability v and the number of trials defined by the number of 
infectious individuals in that subgroup. 

(ii) Commuting process. At each time step, after all individuals update the con- 
tagion dynamics, the simulation moves to their commuting process. The number of 
susceptible travelers departing from their resident subpopulation x(y) per unit time is 
generated from a binomial distribution with probability a and the number of trials 
S xx (t)(S yy (t)). The number of susceptible commuters returning to their residence x(y) 
per unit time is also generated from a binomial distribution with probability t and the 
number of trails S xy (t)(S yx (t)). The number of infectious individuals leaving from or 
returning to their residence is acquired in the same way. 
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